Marine Protected Area Activity

Some R packages will need to be installed so that we can see exactly what is going in with our data. Run the code chunk below to install the packages necessary for this activity. If the package is already installed, the script will only load the package without re-installing it.

## First specify the packages of interest
packages = c("sf", "robis","RColorBrewer",
             "ggplot2", "dplyr", "leaflet")

## Now load or install&load all
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

The U.S. has established nearly 1,000 Marine Protected Areas to protect important places in our ocean, estuaries, coastal waters, and Great Lakes. Marine Protected Areas strive to combat climate change, conserve cultural heritage, and protect marine organisms.

Marine protected areas allow for ecological connectivity, and serve as protected “stepping stones” for marine organisms. Let’s make a map so that we can see all of these stepping stones off of the Coast of California.

First, we’ll need to load in some data that has information on all of the Marine Protected Areas in the United States. The data is in the form of a “shape file”. We will load this file and call it “DataBase”.

# loading our database by setting the directory where the shape file lives with the file extension (.shp)
DataBase = st_read("C:/Users/nposd/Documents/Outreach/MPAs/mpa_inventory_2014_public_shp.shp")

Reading layer `mpa_inventory_2014_public_shp' from data source `C:\Users\nposd\Documents\Outreach\MPAs\mpa_inventory_2014_public_shp.shp' using driver `ESRI Shapefile'
Simple feature collection with 1623 features and 24 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -15.38614 xmax: 180 ymax: 74.70742
Geodetic CRS:  WGS 84

For this example, let’s look at the Marine Protected areas off the Coast of California. When looking at our database, we see that there is a column titled “State”, and each row contains the state that the Marine Protected Area is in. We will filter our database to pull out all of the Marine Protected Areas that are in the “State” of “California” using the following code chunk. We will call this new database “California”.

California = filter(DataBase, State == "California")

There are 185 Marine Protected Areas in our California database! Let’s make a map to see all of the boundaries of the Marine Protected Areas. The “geometry” column has the information of the boundaries. The geometry will give us the latitude and longitude, which we will use to plot the boundary of each Marine Protected Area on our map.

To set up the map, we will first map a blank map of the state of California. We will use built-in “leaflet” library to make our map. Let’s set it up so that we can see the ocean topography using the Esri Ocean Basemap. We’ll save our map as “CaliforniaMap”. To print the CaliforniaMap, we’ll type in CaliforniaMap.

CaliforniaMap = leaflet() %>% # start an empty map
   addProviderTiles("Esri.OceanBasemap") %>% # add the ocean basemap 
   fitBounds( -124, 41, -116, 32) # zoom in to the california area

CaliforniaMap # print the map

Nice! Now let’s add the boundaries of all of the Marine Protected Areas that are on the California Coast. We will loop through each individual Marine Protected area, and plot the latitudes and longitudes on our California map, that way, we can see ALL of the Marine Protected Areas on our map together. We will use “addPolygons” to add each Marine Protected Area’s boundary. You can pick the color that you would like to make the boundary in the “color” section of “addPolygons”. For this example, I made the color red, but you can change it to whatever color you would like!

# loop through each MPA and plot the boundaries on our California Map
for(i in 1:length(California$Site_Name)) {
Geo = as.data.frame(California[[25]][[i]][[1]][[1]]) #the geometry is in the 25th column
Longitude = Geo$V1 #Longitude is saved as Variable 1
Latitude = Geo$V2 #Latitude is saved as Variable 2
CaliforniaMap = CaliforniaMap  %>% 
       addPolygons(lng=Longitude,lat=Latitude, 
                   color = "red", # I wanted to make the MPAs red, but you can put any colors you want!
                   opacity = .5, 
                   fillOpacity = 0.2)
}

CaliforniaMap # print the map

Wow, California has a LOT of Marine Protected Areas! Let’s analyze some data from Marine Protected Areas in California to see exactly how they have helped in the protection of marine organisms! We will need to connect to an online databases to access information about marine organisms.

Let’s choose a specific Marine Protected Area we want to learn more about.

SiteName = ('Soquel Canyon State Marine Conservation Area')  #Look at the California data frame in your global environment...
# and choose the MPA you want to learn more about
ChosenMPA = filter(California, Site_Name == SiteName) #filter for that specific MPA
#get the MPA polygon ready
ChosenMPA_hull = st_convex_hull(ChosenMPA$geometry)
ChosenMPA_text = st_as_text(ChosenMPA_hull)

Introduce OBIS here.

SpeciesOccurence = occurrence(geometry = ChosenMPA_text)

Retrieved 5000 records of approximately 7322 (68%)
Retrieved 7322 records of approximately 7322 (100%)

We’re interested in the number of unique species that occur in the chosen MPA. Let’s group our data to find out how many of each species was observed.

# Let's sum up the unique species that occur in the chosen MPA
if ("individualCount" %in% colnames(SpeciesOccurence)){ #if individualCount is an available column..
  SpeciesOccurence$individualCount <- suppressWarnings(as.numeric(SpeciesOccurence$individualCount)) #introduces NAs
  SpeciesOccurence$individualCount[is.na(SpeciesOccurence$individualCount)] <- 1 #convert NANs to 1; I'm assuming that if it's listed, there was at least one count
  SpeciesOccurence$Count <- 1 * SpeciesOccurence$individualCount #make a new column for Counts of each species
} else {
  SpeciesOccurence$Count = 1 #if individualCount is not a column, I'm assuming that if it's listed there was at least 1 count
}

It’s time to start visualizing our data!

First let’s look at what time period our samples were from by plotting a histogram of the years the samples were taken.

ggplot(SpeciesOccurence, aes(x=date_year)) + geom_histogram(binwidth = 1)+
  labs(title="Histogram of Observation Year", x="Year", y = "Count")

Now let’s look at what depth our samples were from by plotting a histogram of the maximum depth the samples were taken.

ggplot(SpeciesOccurence, aes(x=maximumDepthInMeters)) + geom_histogram(binwidth = 10)+
  labs(title="Histogram of Sample Depthr", x="Depth (meters)", y = "Count")

We’ll start with visualizing the different phyla that exist within the MPA. A phylum is a level of classification or taxonmic rank that allows scientists to categrozie animals. Mammals are in the chordate phylum while insects, spiders, millipeds, and crustaceans are in the arthropod phylum. Let’s look to see what phyla exist in this MPA.

Let’s look at the top then animals that were observed in the MPA.

animals = filter(SpeciesOccurence, kingdom == "Animalia") #filtering the dataframe to pull out the animals
TenAnimals = data.frame(sort(table(animals$scientificName), decreasing = TRUE)[1:10]) #finding the 10 animals that occur most often in the dataframe
print(TenAnimals)

Now let’s take a look at the most common animal in the MPA and see when it occured.

MostCommon_Animal = filter(SpeciesOccurence, scientificName == TenAnimals$Var1[[1]]) #filter for the most commonly occuring animal

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls
5: In readChar(file, size, TRUE) : truncating string with embedded nuls
6: In readChar(file, size, TRUE) : truncating string with embedded nuls
7: In readChar(file, size, TRUE) : truncating string with embedded nuls
8: In readChar(file, size, TRUE) : truncating string with embedded nuls

#Plot the most commonly occuring animal in the MPA over time to see when it occured and how many of them were observed at that time
ggplot(MostCommon_Animal, aes(x = date_year)) + 
  geom_histogram(color = "darkblue", fill = "lightblue") + 
  labs(title = "Histogram of Most Commonly Occuring Animal in the MPA", x = "Year", y = "# of Kelp Bass") # make a histogram showing the most commonly occuring animal